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Abstract — A  modification  of  the  statistical  technique  known 
as  Expectation  Maximization  can  he  applied  for  detecting  and 
characterizing  spectral  activity.  The  novel  technique  adapts  the 
Expectation  Maximization  to  process  histogram  data  using  a  sig¬ 
nal  model.  The  method  provides  signal  clustering  and  parameter 
estimation. 

1.  Introduction 

Conventional  signal  processing  methods  can  be  used  to  gain 
information  of  the  spectral  activity.  Spectral  survey  methods 
can  make  use  filtering,  downconversion,  frequency  and  phase 
locking  techniques  [5],  simple  statistics,  and  neural  networks 
[5].  The  methods  can  exploit  unique  features  or  use  specialized 
processing  techniques  [4]. 

This  paper  is  going  to  treat  spectral  survey  as  a  pattern 
recognition  problem.  The  proposed  method  provides  clustering 
and  parameter  estimation  capabilities.  The  development  of 
the  Spectral  Mixture  Models  begins  with  a  reinterpretation  of 
the  statistical  parameters  in  terms  of  the  telecommunications 
quantities.  The  similarities  between  the  statistics  and  telecom¬ 
munications  are  presented  in  Table  I. 

TABLE  I 

Comparison  between  Statistics  and  Telecommunications 


Statistics 

Telecommunications 

Histogram 

Spectrum 

Clusters 

Signals 

Mean  /r. 

Center  Frequency  / 

Standard  Deviation  cr 

Bandwidth  b 

Mixture  probability  a 

Amplitude-bandwidth  product  a. 

Statistical  sample 

- 

- 

Signal  Sample 

The  frequency  spectrum  provides  multimodal  histogram 
information.  We  wonder  if  a  clustering  algorithm  such  as  the 
Gaussian  Mixture  Models  can  be  use  for  gaining  information 
about  the  activity  in  the  spectrum.  The  idea  sounds  interesting, 
but  it  has  a  problem.  There  is  no  equivalent  concept  of  a  sta¬ 
tistical  sample  in  the  telecommunications  domain.  Fortunately, 
it  will  be  shown  that  it  is  possible  to  adapt  the  algorithm  to 
process  histogram  data.  The  process  will  be  discussed  in  this 
paper. 

We  will  find  out  later  that  there  is  a  second  problem.  The 
Gaussian  model  is  not  adequate  for  representing  telecommuni¬ 
cations  signals.  A  more  suitable  model  for  digital  modulation 
signals  will  be  provided. 


II.  Expectation  Maximization 

The  present  discussion  assumes  that  the  reader  is  famil¬ 
iarized  with  the  derivation  of  the  Gaussian  Mixture  Models 
Algorithm  which  is  an  implementation  of  the  Expectation 
Maximization  using  a  Gaussian  distribution. 

The  expectation  of  the  log-likelihood  Q{0)  is  defined  in 
terms  of  the  conditional  expectation  the  known  data 

vector  y'  ,  the  missing  data  z'  which  identifies  the  cluster 
or  mixture  element,  the  parameter  vector  9  and  the  joint 
probability  of  the  data  and  the  mixture  element  p{y',  z';  6). 

Q{0)  =  E^,^^^,{p{y',z'-,e)}  (1) 

The  joint  probability  can  be  separated  in  two  terms:  the 
model  p(y'\z';0)  and  the  mixture  probability  P{z'). 

p{y',z']e)=p{y'\z'\e)P{z']e)  (2) 

The  parameter  vector  9  contains  the  following  parameters: 
the  mixture  probability  a  =  P{z';  9),  the  mean  vector  p  and 
the  covariance  matrix  S.  Maximizing  the  expectation  with 
respect  to  these  parameters  provides  optimal  density  functions 
that  characterize  the  sample  space. 

9=[p,E,a]  (3) 

In  our  spectral  survey  problem,  the  parameter  vector  con¬ 
sists  of  three  quantities:  the  center  frequency  /,  the  bandwidth 
b  and  the  amplitude-bandwidth  product  a. 

^=lf,b,a]  (4) 

The  method  under  development  can  applied  to  multidimen¬ 
sional  cases.  An  interesting  application  the  application  to  the 
joint  time-frequency  domain.  However,  this  discussion  is  out 
of  the  scope  of  our  present  paper.  For  simplicity,  the  problem  is 
restricted  to  the  one  dimensional  case  which  for  the  Gaussian 
Mixture  Models  consists  at  this  moment  of  fictitious  statistical 
samples. 

y'  =  y'  (scalar)  (5) 
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III.  Histogram  Data 

The  Gaussian  Mixture  Models  will  be  adapted  to  work  with 
histogram  data.  The  approach  will  consider  two  sets.  The  first 
set  Y'  describes  J'  statistical  events  that  represent  our  known 
data.  This  is  our  input  data  y'  =  y'^. 

Y'  =  =  (6) 

A  second  set  Y  is  defined  as  the  histogram  of  the  set  Y' . 
Each  partition  j  of  the  sample  space  is  identified  with  a  bin 
number  yj  and  the  histogram  count  Wj.  This  set  is  what  we 
would  like  our  input  to  be  y  =  yj. 

Y  =  {{y„w,):j  =  l,...,J}  (7) 


resulting  in  an  expression  that  is  similar  to  the  original 
Bayesian  estimate. 


n-\-Wj  —1 


p{yj,Zj  =  k;0)=  —  p{y'^  C  yj\z'^  =  k)P{zj  =  k) 

^  m—n 

(13) 


K 


k=l 


pivj'^  ^)  =  Y1  ^)Pizj  =  k)  (14) 

p{yj,Zj  =  k;0) 


p{zj  =  k\yj]0)  = 


PiVY^) 


(15) 


We  could  assume  that  the  elements  of  y'^  are  sorted  in 
ascending  order.  In  such  case,  the  elements  of  a  bin  yj  can  be 
expressed  by: 

Vj  {i/n5  ■*■5  l}  (^) 

Associated  with  each  bin  yj,  there  is  also  the  missing  data 
Zj  which  contains  information  about  the  mixture  element  or 
cluster  index.  (There  is  also  a  corresponding  variable  for 

each  sample  y^.)  The  cluster  element  is  identified  by  an 
integer  number  from  1  to  K,  where  K  is  the  number  of 
clusters. 

l<Zj<K  (9) 

We  will  assume  that  every  sample  point  within  a  histogram 
bin  shares  the  same  missing  variable.  In  average,  the  statistical 
samples  in  a  certain  vicinity  should  belong  to  the  same  cluster. 
We  can  argue  that  if  the  bin  size  is  small  enough,  this 
assumption  will  be  true. 

P{z'J  =  =  ...  =  =  P{z,  =  k)  (10) 

An  alternative  and  perhaps  less  confusing  notation  for  the 
mixture  probability  is  given  by  equation  11.  This  probability 
can  be  represented  as  a  matrix  with  indexes  for  the  bin  and 
the  cluster  numbers. 

P{j,k)  =  P{zj  =  k)  (11) 

The  model  must  also  be  expressed  in  terms  of  the  histogram 
variables.  This  step  requires  to  define  a  new  conditional 
probability  for  each  bin  p{yj\zj;0)  as  an  average  of  the 
probabilities  of  all  the  statistical  samples  within  a  histogram 
bin. 

n-\-Wj  —  1 

p{yj\zj  =  k;0)  =  ^  Y  P(y'm  C  ypzln  =  k;  0)  (12) 

^  ni—n 

Now,  we  can  proceed  with  the  calculation  of  the  joint 
probability  density  and  the  a  posteriori  probability  needed 
by  the  Expectation  Maximization.  The  term  Wj  is  canceled 


The  histogram  data  will  be  incorporated  into  the  conditional 
expectation  of  the  log-likelihood.  Each  probability  of  the 
histogram  bin  data  is  assumed  to  be  a  independent  event. 

The  variable  Wj  indicates  the  number  of  times  that  the  joint 
probability  shown  in  equation  13  appears  in  the  likelihood 
expression.  This  is  equivalent  to  raise  each  probability  term  to 
the  Wj  power.  The  term  P{zj;  0)™  will  be  treated  as  a  single 
term  for  convenience. 

a  =  P{zf,0)  =  P{zf,0r  (16) 

The  variable  Wj  was  originally  defined  as  an  integer,  but  by 
looking  at  equation  17,  we  notice  that  the  range  of  w  can  be 
extended  to  any  a  real  positive  number.  This  quantity  will  be 
used  to  represent  the  amplitude  of  the  spectrum. 

Q{0)  =  OrP{z;  0)}  (17) 

The  maximization  process  requires  the  calculation  of  the 
derivatives  of  with  respect  to  each  parameter:  /,  b  and  a  and 
finding  the  roots.  The  roots  are  the  new  updated  parameters. 
Einding  the  roots  may  require  implementing  methods  such  as 
the  Newton’s  algorithm. 

VeQ(0)  =  O  (18) 

The  parameter  updates  (M-step)  are  used  to  calculate  the  a 
priori  probabilities  (E-step)  using  the  Bayesian  formulas  that 
have  been  derived.  The  process  is  repeated  until  the  process 
converges  to  some  local  maximum. 

0  «  0max  (19) 

The  parameter  vector  was  relabeled  to  use  the  spectral 
survey  parameters.  After  running  several  trials,  the  new  algo¬ 
rithm  was  unable  to  converge  to  meaningful  parameters.  The 
algorithm  processed  a  mixture  of  MPSK  and  MESK  signals. 
The  equations  and  algorithm  appeared  to  be  implemented 
properly.  Using  a  Gaussian  model  was  not  an  good  choice. 
Adjacent  clusters  interfere  with  each  other  due  to  the  wide 
transition  bandwidth.  The  model  also  lacks  of  a  flat  passband. 
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IV.  Clipped  Gaussian  Model 


Signal  Mixtures 
(19  Iterations) 


Our  telecommunication  signal  requires  a  bandwidth  efficient 
model.  The  model  must  incorporate  a  profile  that  approximates 
an  ideal  hlter  response  in  the  frequency  domain.  Our  hrst 
choice  would  be  to  use  a  Butterworth  hlter  model.  Unfor¬ 
tunately,  the  use  of  this  model  results  in  complex  analytical 
equations  during  the  maximization  process.  Our  desire  is 
to  hnd  bandwidth  efficient  models  that  are  computationally 
efficient. 

A  distribution  known  as  Clipped  Gaussian  will  be  used  for 
our  model  [1].  The  model  provides  a  hat  pass  band  and  a 
steep  transition  bandwidth.  The  model  looks  appealing  because 
the  log-likelihood  simplihes  the  product  of  the  exponential 
function.  The  implementation  of  the  algorithm  using  this 
model  produced  promising  results. 

Where,  N  is  an  even  positive  integer  and  r(x)  is  the 
incomplete  Gamma  function.  Selecting  b  =  2'/2cj  and  N  =  2 
results  in  the  Gaussian  distribution,  but  this  is  not  a  good 
choice.  Our  signal  model  was  built  using  N  =  8. 

V.  Overestimation 

The  number  signals  in  the  spectrum  is  assumed  to  be 
unknown.  We  hope  that  initializing  the  algorithm  with  many 
mixtures  will  result  in  overestimation.  This  will  allow  us  to 
capture  more  details  in  the  frequency  spectrum. 

The  overestimation  in  the  Spectral  Mixture  Models  gener¬ 
ates  suboptimal  solutions  which  will  be  referred  as  subbands. 
The  hnal  solutions  will  be  referred  as  bands  or  signals.  The 
algorithm  will  require  an  additional  step  for  recombining 
adjacent  subbands  that  belong  to  the  same  band.  The  recom¬ 
bination  can  make  use  of  the  parameters  to  determine  the 
proximity  with  two  adjacent  subbands.  Recombination  should 
be  done  when  the  algorithm  has  converged  to  some  suboptimal 
solution. 

VI.  Controlling  Convergence 

The  Clipped  Gaussian  model  has  a  singularity  at  5  =  0. 
To  avoid  convergence  problems,  we  can  execute  a  rule  imme¬ 
diately  after  updating  the  parameters.  The  rule  will  limit  the 
growth  of  b.  The  rule  can  be  implemented  with  a  simple  if 
statement. 

if{b  <  bmin)  then  bj  =  bmm  endif 

Other  useful  way  to  control  the  convergence  of  the  algo¬ 
rithm  is  by  reserving  one  mixture  element  for  noise  floor 
characterization.  The  bandwidth  of  one  mixture  is  forced  to 
be  wide  enough  so  it  is  force  to  converge  to  the  receiver  noise 
floor. 


Fig.  1.  Convergence  of  40  spectral  mixtures  after  19  iterations.  Signals  from 
left  to  right:  QPSK,  tone,  QPSK,  FSK-2  and  noise  floor. 

VII.  Simulations 

Tables  II  to  VII  contain  the  results  of  several  Monte  Carlo 
simulations.  The  simulations  were  conducted  using  combina¬ 
tion  of  center  frequencies:  2,  4  and  6  kHz;  bandwidth:  1  kHz; 
and  various  signal-to-noise  ratios  (SNR):  40,  30,  20,  10,  5  and 
2  dB.  The  algorithm  was  initialized  with  uniformly  spaced 
mixtures.  The  initial  bandwidths  were  uniform  except  for  the 
bandwidth  of  one  mixture  that  is  forced  to  converge  to  the 
noise  floor  of  the  receiver.  The  algorithm  was  stopped  after 
25  iterations.  The  signal  under  test  is  a  QPSK  type  sampled 
at  16  kHz.  The  message  is  a  random  sequence  that  contains 
between  400  and  500  symbols.  The  true  number  of  cycle  is  a 
multiple  of  the  total  number  of  channels  in  the  analysis  hlter 
used  to  estimate  the  frequency  spectrum.  A  total  of  256  hlter 
channels  were  used  in  all  the  simulations. 

VHI.  Conclusion 

Expectation  Maximization  can  be  use  for  spectral  survey  by 
adapting  the  algorithm  in  two  ways.  First,  the  algorithm  needs 
to  process  histogram-like  data  such  as  the  frequency  spectrum 
or  spectrogram.  Second,  the  model  needs  to  be  representative 
of  a  signal.  A  display  of  the  spectral  mixtures  is  shown  in 
Figure  1. 

The  concept  can  be  rehned  to  produce  better  estimates.  Bet¬ 
ter  combination  rules  can  be  implemented.  Initial  conditions 
can  be  modihed  to  acelerate  the  convergence.  It  is  important 
to  notice  that  the  parameter  b  was  not  meant  to  be  a  3dB 
bandwidth  estimate.  A  better  bandwidth  estimates  could  be 
produced  by  mapping  b  to  the  real  3dB  bandwidth. 

Potential  applications  of  this  method  include  cognitive  ra¬ 
dios,  broadband  receivers  and  measurement  equipment. 
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TABLE  II 

Center  frequency  estimates  for  various  SNR  levels.  True 
PARAMETERS  CENTER  FREQUENCY  =  2  KHZ,  BANDWIDTH  =  1  KHZ. 


TABLE  VI 

Bandwidth  estimates  for  various  SNR  levels.  True  parameters 
Center  Frequency  =  4  kHz,  Bandwidth  =  1  kHz. 


Fc  =  2kHz 

Center  Frequency  Estimate  (Hz) 

SNR 

Mean 

Bias 

Std.  Dev. 

2 

2005.48 

-5.48 

437.81 

5 

2057.59 

-57.59 

339.73 

10 

2010.82 

-10.82 

278.32 

20 

2015.63 

-15.63 

166.07 

30 

2024.53 

-24.53 

138.25 

40 

2022.31 

-22.31 

135.74 

Fc  =  AkHz 

Bandwidth  Estimate  (Hz) 

SNR 

Mean 

Bias 

Std.  Dev. 

2 

1255.90 

-255.90 

822.41 

5 

1332.06 

-332.06 

677.71 

10 

1454.18 

-454.18 

450.95 

20 

1652 

-652 

273.95 

30 

1709.18 

-709.18 

233.25 

40 

1738.41 

-738.41 

211.03 

TABLE  III  TABLE  VII 

Center  frequency  estimates  for  various  SNR  levels.  True  Bandwidth  estimates  for  various  SNR  levels.  True  parameters 
PARAMETERS  CENTER  FREQUENCY  =  4  KHZ,  BANDWIDTH  =  1  KHZ.  CENTER  FREQUENCY  =  6  KHZ,  BANDWIDTH  =  1  KHZ. 


Fc  =  AkHz 

Center  Frequency  Estimate  (Hz) 

SNR 

Mean 

Bias 

Std.  Dev. 

2 

4303.15 

-303.15 

70.04 

5 

4026.37 

-26.37 

144.73 

10 

4020.42 

-15.80 

52.56 

20 

4015.80 

-15.65 

8.15 

30 

4015.65 

-15.70 

4.80 

40 

4015.70 

-15.70 

2.91 

Fc  =  6kHz 

Bandwidth  Estimate  (Hz)  | 

SNR 

Mean 

Bias 

Std.  Dev. 

2 

1100.72 

-100.72 

671.86 

5 

1167.12 

-167.12 

593.57 

10 

1380.17 

-380.17 

488.73 

20 

1649.5 

-649.5 

319.96 

30 

1752.25 

-752.25 

248.30 

40 

1764.62 

-764.62 

239.42 

TABLE  IV 

Center  frequency  estimates  for  various  SNR  levels.  True 
PARAMETERS  CENTER  FREQUENCY  =  6  KHZ,  BANDWIDTH  =  1  KHZ. 


Fc  =  6kHz 

Center  Frequency  Estimate  (Hz) 

SNR 

Mean 

Bias 

Std.  Dev. 

2 

5339.43 

660.56 

144.27 

5 

5957.93 

42.06 

42.23 

10 

6011.36 

-14.84 

5.40 

20 

6014.84 

-15.46 

10.32 

30 

6015.46 

-15.68 

4.74 

40 

6015.68 

-15.68 

3.55 
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TABLE  V 

Bandwidth  estimates  for  various  SNR  levels.  True  parameters 
Center  Frequency  =  2  kHz,  Bandwidth  =  1  kHz. 


Fc  =  2kHz 

Bandwidth  Estimate  (Hz) 

SNR 

Mean 

Bias 

Std.  Dev. 

2 

1063.18 

-63.18 

643.80 

5 

1249.49 

-249.49 

603.34 

10 

1409.90 

-409.90 

472.68 

20 

1670.68 

-670.68 

303.80 

30 

1747.25 

-747.25 

255.34 

40 

1796.64 

-760.98 

250.65 
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